\section{Growing Large Dust Grains}
\label{appendix:graph}

In this chapter, we seek a realistic but computationally efficient means
of growing dust grains of arbitrarily large size.  It is possible to use
network calculations that follow all species up to those with $\sim 10^4$
atoms, but following larger species requires a different approach.

To treat this problem, we follow the time evolution of a network
with branchings \cite{wang09}.  In such a case, we analyze the
abundance $Y_i(t + \Delta t)$ of a species $i$ in a linear network
at time $t + \Delta t$ in terms of branchings of a directed graph.  The
vertices of the graph are the species in the network.  The directed arcs
in the graph are the reactions between the species.  Since the network
solutions at $t + \Delta t$ are in terms of the abundances at $t$,
an arc from vertex $i$ to vertex $j$ gives the contribution of vertex
$j$ at time $t$ to the abundance of $i$ at time $t + \Delta$.  As a
consequence of this, and because of the rules of branching weights,
the weight of the arc from $i$ to $j$ is $\ln\left( \lambda_{ji} \Delta t
\right)$, where $\lambda_{ji}$ is the rate of the reaction from $j$ to
$i$.  Note that direction of the arc is opposite that of the reaction in
time.

A branching $B$ in a graph is a set of arcs in the graph such that there
are no cycles in the graph and no vertex has indegree (that is, the
number of arcs entering a vertex) larger than one.  The weight $w(B)$ of a
branching $B$ is the sum of the weights of all the arcs in the branching.

With these preliminaries in mind, the abundance of species $i$ in the
network at $t + \Delta t$ is
\begin{equation}
Y_i(t + \Delta t) = \frac{\sum_{B:(i)} \exp\{w(B)\} Y_B^{(i)}(t)}{\sum_{B}
\exp\{w{B}\}}
\label{eq:xi}
\end{equation}
The symbol $B:(i)$ indicates a branching rooted at $i$, that is, a branching
such that vertex $i$ has no arcs entering it.  The symbol $Y_B^{(i)}(t)$
is defined as
\begin{equation}
Y_B^{(i)}(t) = \sum_{j:\{(i)\to j\}_B} Y_j(t)
\label{eq:xbi}
\end{equation}
where the symbol $j:\{(i)\to j\}_B$ means any vertex $j$ for which there is
a path from root $i$ to $j$ in branching $B$.

We now turn to the problem of forming large dust grains.  We consider
a network that may be split into two subnetworks--one has complicated
interactions among the subnetwork members and the other has a simple
capture chain.  Fig. \ref{fig:net} shows a simple six species example.
Species 1 through 3 interact via capture
and destruction reactions.  Species 4 through 6 only interact via capture
reactions.  We may thus split up the network into two subnetworks, $C_1$ and
$C_2$, and the parent digraph for branchings is then that in
Fig. \ref{fig:2-cluster}.

\begin{figure}
\centering
\includegraphics[height=1in]{figures/net}
\caption{A simple 6 species network.}
\label{fig:net}
\end{figure}

\begin{figure}
\centering
\includegraphics[height=1in]{figures/2-cluster}
\caption{A simple 6 species network.}
\label{fig:2-cluster}
\end{figure}

We now consider the general case.  We imagine the $n$ species in $C_1$ to
be labeled 1 through $n$ and the species in $C_2$ to be labeled $n+1$ to
$N$.  We also consider that $C_1$ and $C_2$ are linked via a capture
from $n$ to $n+1$.
If we focus on a subnetwork consisting of $C_2$ plus species $n$, the possible
exponential sum of branchings is
\begin{equation}
\sum_{B_2} \exp\{w({B_2})\} =
\prod_{i = n}^N \left( 1 + \lambda_i \Delta t\right)
\label{eq:wb2}
\end{equation}
because each reaction in the subnetwork will either contribute an arc (weight
$\ln \left(\lambda_j \Delta t\right)$) or not (weight zero) to a given
branching.
Consider now the full network.  We build up branchings from a sub-branching
within $C_1$ and within $C_2$ and links between them.
If a branching contains the arc from
$C_2$ to $C_1$, that is, from $n+1$ to $n$,
the sub-branch for the subnetwork $C_1$ must be rooted at $n$ because
there can only be one arc entering $n$.  On the other hand, if the
branching does not include the arc from $n+1$ to $n$, the sub-branching within
$C_1$ can be rooted at any of the vertices in $C_1$.
The exponential sum of all branchings will then be
\[
\sum_B \exp\{w(B)\} = \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
\prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
+
\sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
\prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
\]
\begin{equation}
= \left( \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
  \right)
  \prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
\label{eq:btotal}
\end{equation}
where $B_1:[(k)]_{C_1}$
indicates the branching $B_1$ is rooted at $k \in C_1$ and is entirely
contained in $C_1$.

We begin by considering the abundance of species $n$.  The abundance
of species $n$ at $t + \Delta t$ is given by Eq. (\ref{eq:xi}).  Since
the relevant branchings must be rooted at $n$, they cannot contain any
arcs entering vertex $n$; thus, we find
\begin{equation}
Y_n(t + \Delta t) =
\frac{
  \sum_{B_1:[(3)]_{C_1}} \exp\{w(B_1)\}
  \prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t \right)
  Y_{B_1}^{(n)}
}{
\left( \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
  \right)
  \prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
}
\end{equation}
This reduces to
\begin{equation}
Y_n(t + \Delta t) =
\frac{
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\}
  Y_{B_1}^{(n)}
}{
 \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
}
\label{eq:xn}
\end{equation}

We now consider the abundance of species $j > n$.
Since the relevant branchings must be rooted at $j$, they cannot contain the
arc entering vertex $j$.
The other question we address about the relevant
branching is whether it includes the from $n+1$ to $n$.
Branchings that do not include that arc will include the factor
$\sum_{B_1:[(k\in C_1)]_{C_1}} \exp\{w(B_1)\}$.  Branchings that do
include that arc will include the factor
$\sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\}\lambda_n \Delta t$.
The contribution to the numerator in Eq. (\ref{eq:xi}) will be
\begin{equation}
N_j =
\sum_{l = n+1}^j Y_l(t)
\left( \prod_{m = n+1}^{l-1} \left( 1 + \lambda_m \Delta t \right) \right)
\left( \prod_{p = l}^{j-1} \lambda_p \Delta t \right)
\left( \prod_{q = j+1}^N \left( 1 + \lambda_q \Delta t \right) \right)
\label{eq:num_part1}
\end{equation}
because the branching must have a path from $l$ to $j$ for $Y_l(t)$ to
contribute to $Y_j(t + \Delta t)$.
The result is that
\[
Y_j(t + \Delta t) =
\frac{
\left(
  \sum_{B_1:[(k\in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\}\lambda_n \Delta t
\right) N_j
}
{
\left( \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
  \right)
  \prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
}
\]
\begin{equation}
+
\frac{
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\}\lambda_n \Delta t Y_{B_1}^{(n)}
  \left( \prod_{p = n+1}^{j-1} \lambda_p \Delta t \right)
}
{
\left( \sum_{B_1:[(k \in C_1)]_{C_1}} \exp\{w(B_1)\}
  +
  \sum_{B_1:[(n)]_{C_1}} \exp\{w(B_1)\} \lambda_n \Delta t
  \right)
  \prod_{i = n+1}^N \left( 1 + \lambda_i \Delta t\right)
}
\label{eq:xj1}
\end{equation}
This reduces to
\begin{equation}
Y_j(t + \Delta t) =
\sum_{l = n+1}^j \frac{Y_l(t)}{1 + \lambda_j \Delta t}
\left( \prod_{p = l}^{j-1}
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}\right)
+
Y_n(t + \Delta t)
\left(
  \frac{\lambda_n \Delta t}{1 + \lambda_j \Delta t}
\right)
\prod_{p = n+1}^{j-1} \left(
   \frac{
      \lambda_p \Delta t
   }{
      1 + \lambda_p \Delta t
   }\right)
\label{eq:xj2}
\end{equation}
When $\Delta t \to 0$, we see $Y_j(t + \Delta t) \to Y_j(t)$ while when
$\Delta t \to \infty$, we find $Y_j(t + \Delta t) \to
\frac{\lambda_n}{\lambda_j}Y_n(t + \Delta t)$.

Eq. (\ref{eq:xj2}) gives the abundance of a species $j > n$ at $t + \Delta t$
in terms of species $Y_n(t + \Delta t)$ and species with index $k < j$ at
$t$.  We now consider allowing only certain species with $j > n$ to have
non-zero abundance.  We label these $j_1$, $j_2$, and so forth.  With
our assumption, we insist that species with index $k$ such that $n < k < j_1$
have zero abundance.  We thus find
\begin{equation}
Y_{j_1}(t + \Delta t) = \frac{Y_{j_1}(t)}{1 + \lambda_{j_1}\Delta t}
+ Y_n(t + \Delta t)\left(\frac{\lambda_n\Delta t}{1 + \lambda_{j_1}\Delta t}
\right)\prod_{p=n+1}^{j_1-1}\left(\frac{\lambda_p \Delta t}
{1 + \lambda_p \Delta t}\right)
\label{eq:xj_1}
\end{equation}
Similary, species with index $k$ such that $j_1 < k < j_2$ have zero abundance;
thus,
\[
Y_{j_2}(t + \Delta t) =
\frac{Y_{j_1}(t)}{1 + \lambda_{j_2}\Delta t}
\prod_{p=j_1}^{j_2-1}\left(\frac{\lambda_p \Delta t}
{1 + \lambda_p \Delta t}\right)
+
\frac{Y_{j_2}(t)}{1 + \lambda_{j_2}\Delta t}
\]
\begin{equation}
+
Y_n(t + \Delta t)\left(\frac{\lambda_n\Delta t}{1 + \lambda_{j_2}\Delta t}
\right)\prod_{p=n+1}^{j_2-1}\left(\frac{\lambda_p \Delta t}
{1 + \lambda_p \Delta t}\right)
\label{eq:xj_2}
\end{equation}
This can be written
\begin{equation}
Y_{j_2}(t + \Delta t) =
\frac{1 + \lambda_{j_1}\Delta t}{1 + \lambda_{j_2}\Delta t}
\prod_{p = j_1}^{j_2 - 1}\left(
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
\right) Y_{j_1}(t + \Delta t)
+
\frac{Y_{j_2}(t)}{1 + \lambda_{j_2}\Delta t}
\label{eq:xj2b}
\end{equation}
which, in turn, becomes
\begin{equation}
Y_{j_2}(t + \Delta t) =
\frac{\lambda_{j_1}\Delta t}{1 + \lambda_{j_2}\Delta t}
\prod_{p = j_1 + 1}^{j_2 - 1}\left(
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
\right) Y_{j_1}(t + \Delta t)
+
\frac{Y_{j_2}(t)}{1 + \lambda_{j_2}\Delta t}
\label{eq:xj3b}
\end{equation}
A similar analysis holds for subsequent species such that
\begin{equation}
Y_{j_m}(t + \Delta t) =
\frac{\lambda_{j_{m-1}}\Delta t}{1 + \lambda_{j_m}\Delta t}
\prod_{p = j_{m-1} + 1}^{j_m - 1}\left(
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
\right) Y_{j_{m-1}}(t + \Delta t)
+
\frac{Y_{j_m}(t)}{1 + \lambda_{j_m}\Delta t}
\label{eq:xjm}
\end{equation}
Note that it is important to distinguish between the symbols
$j_{m-1}$, which is the $(m-1)^{th}$ non-zero species in our approximation,
and $j_m - 1$, which is the species with index one less than $j_m$.
We are in fact interested in an effective rate for capture from $j_{m-1}$
to $j_m$, so we compute
\[
\frac{
Y_{j_m}(t + \Delta t) - Y_{j_m}(t)
}
{
\Delta t
}
=
\]
\begin{equation}
\frac{\lambda_{j_{m-1}}}{1 + \lambda_{j_m}\Delta t}
\prod_{p = j_{m-1} + 1}^{j_m - 1}\left(
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
\right) Y_{j_{m-1}}(t + \Delta t)
-
\frac{\lambda_{j_m}}{1 + \lambda_{j_m}\Delta t} Y_{j_m}(t + \Delta t)
\label{eq:dxjm}
\end{equation}
From this we infer an effective capture rate from $j_{m-1}$ to $j_m$ as
\begin{equation}
\lambda_{j_{m-1},j_m}^{eff} =
\frac{\lambda_{j_{m-1}}}{1 + \lambda_{j_m}\Delta t}
\prod_{p = j_{m-1} + 1}^{j_m - 1}\left(
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
\right)
\label{eq:lambda_jm_eff}
\end{equation}

Eq. (\ref{eq:dxjm}) is still incomplete in our treatment because the
capture term out allows for a single capture.  This is inconsistent with
our assumption that species with index $k$ such that $j_m < k < j_{m+1}$
have zero abundance.  We thus modify Eq. (\ref{eq:dxjm}) to read
\begin{equation}
\frac{
Y_{j_m}(t + \Delta t) - Y_{j_m}(t)
}
{
\Delta t
}
=
\lambda_{j_{m-1},j_m}^{eff} Y_{j_{m-1}}(t + \Delta t)
-
\lambda_{j_m,j_{m+1}}^{eff} Y_{j_m}(t + \Delta t)
\label{eq:dxjm2}
\end{equation}
This reduces the problem to a simple capture chain with effective capture
rates given by Eq. (\ref{eq:lambda_jm_eff}).

Since the network is nonlinear, smaller species are consumed to create larger
species.  If we imagine capture is purely of single atoms and that
larger species are made of a single type of atomic species $a$, the contribution
to the change
in the abundance of single atoms when species $j_{m-1}$ captures to $j_m$
will be
\begin{equation}
\frac{
Y_{a}(t + \Delta t) - Y_{a}(t)
}
{
\Delta t
}
=
-
\left( N_{j_m} - N_{j_{m-1}} \right)
\left(
\lambda_{j_{m-1},j_m}^{eff} Y_{j_{m-1}}(t + \Delta t)
-
\lambda_{j_m,j_{m+1}}^{eff} Y_{j_m}(t + \Delta t)
\right)
\label{eq:dxa}
\end{equation}
where $N_{j_{m-1}}$ is the number of atoms $a$ that make up species
$j_{m-1}$ and $N_{j_m}$ is the number of atoms $a$ that make up species
$j_m$.

To solve the network equations by the Newton-Raphson method, we define
\begin{equation}
f_{j_m} = 
\frac{
Y_{j_m}(t + \Delta t) - Y_{j_m}(t)
}
{
\Delta t
}
-
\lambda_{j_{m-1},j_m}^{eff} Y_{j_{m-1}}(t + \Delta t)
+
\lambda_{j_m,j_{m+1}}^{eff} Y_{j_m}(t + \Delta t)
\label{eq:f_jm}
\end{equation}
These (and the other coupled equations) are iterated on updated solutions
at $t + \Delta t$ until convergence (all $f's = 0$).  We compute the
necessary matrix elements as
\begin{equation}
A_{j_m,j_m} = \frac{\partial f_{j_m}}{\partial Y_{j_m}(t + \Delta t)}
= \frac{1}{\Delta t} + \lambda_{j_m,j_{m+1}}^{eff}
\label{eq:amm}
\end{equation}
and
\begin{equation}
A_{j_m,j_{m-1}} = \frac{\partial f_{j_m}}{\partial Y_{j_{m-1}}(t + \Delta t)}
= - \lambda_{j_m,j_{m-1}}^{eff}
\label{eq:ammm}
\end{equation}
We also find the contribution to the right-hand-side vector
$b_{j_m} = -f_{j_m}$.
As for the single atoms $a$, the contribution to the function $f_a$
is
\begin{equation}
f_a =
\frac{
Y_{a}(t + \Delta t) - Y_{a}(t)
}
{
\Delta t
}
+
\left( N_{j_m} - N_{j_{m-1}} \right)
\left(
\lambda_{j_{m-1},j_m}^{eff} Y_{j_{m-1}}(t + \Delta t)
-
\lambda_{j_m,j_{m+1}}^{eff} Y_{j_m}(t + \Delta t)
\right)
\label{eq:fa}
\end{equation}
The contribution to the matrix element $A_{a,j_m}$ is
\begin{equation}
A_{a,j_m} = \frac{\partial f_a}{\partial Y_{j_m}(t + \Delta t)}
= - \left( N_{j_m} - N_{j_{m-1}} \right)
\lambda_{j_m,j_{m+1}}^{eff}
\label{eq:aam}
\end{equation}
and to $A_{a,j_{m-1}}$ is
\begin{equation}
A_{a,j_{m-1}} = \frac{\partial f_a}{\partial Y_{j_{m-1}}(t + \Delta t)}
= \left( N_{j_m} - N_{j_{m-1}} \right)
\lambda_{j_{m-1},j_m}^{eff}
\label{eq:aamm}
\end{equation}
The contribution to $b_a$ is $-f_a$.  Once all matrix and vector elements
are available, the matrix equation $A \delta Y = b$ is solved, the
vector $Y(t + \Delta t)$ is updated as $Y + \delta Y$ and the procedure is
repeated until $\delta X = 0$.

To implement this procedure, we must compute the effective rate in
Eq. (\ref{eq:lambda_jm_eff}).  We first note we will take an individual reaction
rate $\lambda_p$ to be given by $\lambda_p = N_p^{2/3} \lambda_0$.
$N_p^{2/3}$ is proportional to the surface area of species $p$.  $\lambda_0$
itself depends on the abundance of free atoms and other factors like the
sticking probability, which we assume to be the same for all species.
With this in mind, we make the following approximation:
\begin{equation}
\frac{\lambda_p \Delta t}{1 + \lambda_p \Delta t}
= \frac{1}{1 + 1 / \lambda_p \Delta t} \approx
\exp\left(-\frac{1}{\lambda_p \Delta t}\right)
\label{eq:exp_approx}
\end{equation}
With this approximation, we may compute the effective rate in
Eq. (\ref{eq:lambda_jm_eff}) as
\begin{equation}
\lambda_{j_{m-1},j_m}^{eff} =
\frac{\lambda_0 N_{j_{m-1}}^{2/3}}{1 + \lambda_0 N_{j_m}^{2/3}\Delta t}
\exp\left(-\frac{1}{\lambda_0\Delta t}
\sum_{p = j_{m-1}}^{j_m - 1} 
\frac{1}{N_p^{2/3}}\right)
\label{eq:lambda_jm_eff_approx}
\end{equation}
Since the number $N_p$ in the sum in the exponential of
this equation is typically a very large number, we approximate the
sum as an integral:
\begin{equation}
\sum_{p = j_{m-1}}^{j_m - 1} \frac{1}{N_p^{2/3}}
\approx
\int_{N_{j_{m-1}}}^{N_{j_m - 1}} \frac{dN}{N^{2/3}}
=
3 \left(N_{j_m - 1}^{1/3} - N_{j_{m-1}}^{2/3}\right)
\label{eq:sum_approx}
\end{equation}
which finally yields
\begin{equation}
\lambda_{j_{m-1},j_m}^{eff} =
\frac{\lambda_0 N_{j_{m-1}}^{2/3}}{1 + \lambda_0 N_{j_m - 1}^{2/3}\Delta t}
\exp\left(-3\frac{\left(N_{j_m - 1}^{1/3} - N_{j_{m-1}}^{1/3}\right)}
{\lambda_0\Delta t}\right)
\label{eq:lambda_jm_eff_approx_2}
\end{equation}

The individual rate, in our species per atom notation, is
$\lambda_0 = N_A \langle \sigma_0 v \rangle \rho Y_a$.
ADD SOME TEXT HERE.
When $\lambda_0 \Delta t$ is large compared to $N_{j_m}^{1/3}$, the
exponential term does not cut off the effective rate significantly.
When $\lambda_0 \Delta t$ is small compared to $N_{j_m}^{1/3}$ because
either the density or $Y_a$ is small, the exponential will cutoff the
effective rate and limit the build up of heavy species.

